Computer simulations on the sympatric speciation modes for the Midas cichlid 

species complex 



r» ■ 
o ' 

o : 

wo: 
s ■ 

<: 

<n : 
<n . 

E5 1 : 
cm: 

6 : 

x> : 



> 
O 

cn 

od 
o 

o 



X 



K. Luz-Burgoa, 1 S. Moss de Oliveira, 1 and J. S. Sa Martins 1 

1 Instituto de Fisica, Universidade Federal Fluminense, 
Campus da Praia Vermelha, Boa Viagem, Niteroi, 24210-340, RJ, Brazil^ 
(Dated: February 1, 2008) 

Cichlid fishes are one of the best model system for the study of evolution of the species. Inspired 
by them, in this paper we simulated the splitting of a single species into two separate ones via 
random mutations, with both populations living together in sympatry, sharing the same habitat. 
We study the ecological, mating and genetic conditions needed to reproduce the polychromatism and 
polymorphism of three species of the Midas Cichlid species complex. Our results show two scenarios 
for the A. Citrinellus speciation process, one with and the other without disruptive natural selection. 
In the first scenario, the ecological and genetic conditions are sufficient to create two new species, 
while in the second the mating and genetic conditions must be synchronized in order to control the 
velocity of genetic drift. 



I. INTRODUCTION 

The evolution of a single population into two or more 
species without prevention of gene flow through geo- 
graphic segregation is known as sympatric speciation 
[11)13) 9- Cichlid fishes, Amphilophus Zaliosus, are one 
possible example of evolution by sympatric speciation in 
nature [1, Q • The Midas cichlid species complex [|| are 
distributed in the Great Lakes of Nicaragua as well as in 
several crater lakes in the area. Combined, they represent 
by far the largest biomass of any fish species in Nicaragua 
freshwaters. There are substantial morphological differ- 
ences between them, for example, some of them have 
cryptic colouring, grey or brown with dark bars or spots, 
known as Normal Morph, or a conspicuous form, which 
lacks melanophores, resulting in brightly red, orange, or 
yellowish coloured fish, Gold Morph. Cichlid fishes are 
also characterized by a pair of jaws in the pharyngeal 
area in addition to the oral jaws, and this key innovation 
is presumed to be responsible for their great ability to 
colonize new habitats and to exploit successfully a large 
diversity of trophic niches. 

Three different species have already been recognized 
!(| within the Midas cichlid complex: 1.- Amphilophus 
Labiatus, the red devil cichlid, a fleshy-lipped species 
thought to be restricted to the big lakes Lake Nicaragua 
and Lake Managua, 2.- A. Zaliosus, the arrow cichlid, 
an elongated species that is restricted in its distribution 
to one of the crater lakes, Lake Apoyo, and 3.- Amphilo- 
phus citrinellus, the Midas cichlid, a generalist species 
with very widespread distribution. For A. citrinellus, 
a polychromatism has been described Q, with normal 
and gold morphs. Strong assortative mating according 
to colour has been observed as well H, both in the field 
and in captivity, suggesting that sexual selection main- 
tains the colour polymorphism. Two types of pharyngeal 



jaws of A. citrinellus have been described [9|, a papilli- 
form morph with slender pointed teeth and a mollariform 
morph with thicker rounded teeth. These previous stud- 
ies showed that there is a trade-off in performance: the 
mollariform fish are specialized and more efficient at eat- 
ing hard diets such as snails, whereas they are less effi- 
cient in feeding on soft diets than the papilliform morphs 
and vice versa. For A. labiatus, the same polychromatism 
has been described 0] , with normal and gold morphs, and 
only papilliform pharyngeal jaws have been documented. 
All A. zaliosus have cryptic colouring, normal morph, 
and are polymorphic with papilliform and mollariform 
pharyngeal, Fig. [TJ 
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FIG. 1: The scheme shows the polycromatism and polymor- 
phism of the three species of the Midas cichlid species complex 
from the lakes of Nicaragua described. 

According to theoretical models [lj|, sympatric speci- 
ation is driven by disruptive frequency-dependent natu- 
ral selection, caused by competition for diverse resources 
pdl Ojj]. On the other hand, some authors have argued 
that sexual selection can also cause sympatric speciation 
[l3L [3] • Cichlids are renowned for their vast diversity of 
trophic morphologies and often extreme degree of ecologi- 
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cal specialization [9(. However, the sympatric occurrence 
of many sibling species that seem to differ only in colour- 
ing makes it unlikely for ecological specialization to be 
the sole mechanism of speciation in this group. At the 
same time, genetic data (J) show that sympatric specia- 
tion by sexual selection alone is rather unlikely for the 
speciation case of the crater lake Apoyo. We studied the 
ecological, mating and genetic conditions needed to re- 
produce the polychromatism and polymorphism of the 
three species of the Midas Cichlid species complex, Fig. 
[TJ Our study was based on simulations of an individual- 
based model where natural selection caused by compe- 
tition for diverse resources and sexual selection, tuned 
by strength parameters on two quantitative and inde- 
pendent traits, were considered. Our results show two 
scenarios for the sympatric speciation of A. Citrinellus 
species, one with and the other without disruptive natu- 
ral selection. In the first scenario, A. Zaliosus develops 
jaw polymorphism while retaining a single colour morph, 
while in the second A. Labiatus develops polychromatism 
with a single jaw morph. 



II. COMPUTATIONAL MODEL 



The individuals' genomes are represented by three 
pairs of bitstrings, each of them consisting of a computer 
word of 32 bits. The first pair is age-structured and con- 
tains the information of when, between 1 and 32 time 
steps, the individual would die if only genetic diseases 
were considered [l7| . The second pair delineates the abil- 
ity of the individual to survive under a competition for 
the available resources, representing an ecological trait 
such as the types of pharyngeal jaws in the case of cich- 
lid fishes. The third pair represents a trait only related 
to sexual selection, a mating trait, such as the colour 
of cichlid fishes. At the beginning of the simulation, all 
individuals are born with random genomes. When a fe- 
male succeeds in staying alive until reaching a minimum 
reproduction age, A, it looks for a male to mate with and 
generates b offspring every time step before dying, with a 
new choice for a mate being done at each time step. The 
first pair of the offspring's genome is constructed in the 
following way: each one of the first pair of strings of the 
male, for instance, is broken at the same random position 
and the complementary pieces, originated from different 
strings, are joined to form two male gametes. One of the 
gametes is then randomly chosen to be passed on to the 
offspring. After that, one random bad mutation is intro- 
duced into this gamete, and the final result corresponds 
to one string of the new individual. The other string 
is constructed from the first pair of the female's strings 
by the same process, that simulates random crossover, 
recombination and addition of one bad mutation. 



A. Natural selection caused by competition 

In the present model, competition for food is related 
to a phenotype, j, represented by the second, non age- 
structured, pair of bit-strings, which is constructed in 
the same way as the first pair of the individual' genome. 
This phenotypic characteristic is computed by counting 
the number of bit positions where both bits are set to 
1, plus the number of dominant positions (chosen as 16) 
with at least one of the two bits set. It will therefore 
be a number j between and 32, which we will refer to 
as the individual's phenotype. We call Mj the mutation 
probability per locus of this ecological trait. A mutation 
can change the locus either from to 1 or from 1 to 0. 
In order to control the population's size and introduce a 
competition we used the Vcrhulst factor, V(j, t). We con- 
sidered three intra-specific competitions, 16|, depending 
on the individual's phenotype j, each one related to a 
given phenotypic group: 

Vi(j,t), < j < ni; specialist, 
= \ V m (j,i), rii < j < n 2 ; intermediate, 
V 2 (j,t), n 2 < j < 32; specialist., 

where n\ and n-i are two parameters of the simulation. 
At every time step, t, and for each individual with pheno- 
typic characteristic j, a random real number uniformly 
distributed between and 1 is generated; if this num- 
ber is smaller than V(j,t), the individual dies. For the 
specialist groups the competition is given by: 



v m (j,t) = _ , 



(2.2) 



where P\{2) (t) accounts for the population with pheno- 
type j < n\ (j > n-i) at time t, P m (t) accounts for the 
population with phenotype j 6 [711,712], and F(j,t) is a 
resource distribution. Individuals with intermediate phe- 
notypes (P m ) compete among themselves and also with 
a fraction X of each specialist population. The Vcrhulst 
factor for them is: 



V m (j,t) = 



P m (t)+Xx [P^t) + P 2 (t)} 



(2.3) 



Eq. (|2.2p means that specialist individuals (Pi, P 2 ) com- 
pete with those belonging to the same phenotypic group 
and also with the whole intermediate population, but 
there is no competition between specialists of different 
groups because we are assuming they are specialized to 
some extent ([0, ni),(n2, 32]) on particular resources, as is 
the case for papilliform and mollariform pharyngeal jaws 
in the Midas cichlid species complex. In equations 12.21 
and 12.31 the resource distribution used varies according 
to: 



F(j,t) = Cx (1- GOO), with 
G(j) = Zxe-M 2 /« 



(2.4) 
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where C is a carrying capacity, and for all simulations it 
was set to C = 2 x 10 5 . The first case, Z = in Eq. (j2~lj) . 
is used to simulate a scenario without disruptive selection 
on the ecological trait. For this case, all individuals have 
the same carrying capacity, Eq. (|2.2p and Eq. (|2.3p . 
The second case, Z > in Eq. (|2.4j) . is used to simu- 
late disruptive natural selection with a strength Z. For 
this case, all individuals with intermediate phenotypes 
are disadvantaged, with respect to specialist individuals, 
according to a reversed gaussian, G(k). 



B. Sexual selection 

In the simulations, sexual selection is related to an- 
other phenotype, k, and was represented by a new pair of 
non age-structured bit-strings, the mating trait, that also 
obeys the general rules of crossing and recombination. 
This phenotype was computed in the same way as that 
for the ecological trait. It will therefore be a number k 
between and 32, and we call the mutation probabil- 
ity per locus of this mating trait, which can also mutate 
back and forth between and 1. In order to consider as- 
sortative mating in a sympatric environment, we defined 
two phenotypic groups, one composed by the individu- 
als that have k G [0, 16) while individuals of the other 
have k G (16,32]. With some probability, Y G [0.0, 1.0], 
a female with phenotype k will mate with a male of the 
same phenotypic group and with probability 1.0 — Y will 
mate with a male of the other phenotypic group, at each 
time step of its life. For instance, if a female has pheno- 
type k G [0, 16) and a random real number, r, uniformly 
distributed between and 1, is tossed that is smaller 
than Y, it selects its mate among N m males of the same 
phenotypic group k G [0, 16) by picking the one with 
the smallest phenotype value k. If a female has pheno- 
type k G [0, 16) and the random real number, r, is larger 
than Y, it chooses a partner from the other phenotypic 
group, k G (16,32]. A similar rule applies to females 
with k G (16,32], with the proviso that now the female 
picks as mate, among N m males of the same phenotypic 
group, the one with the largest phenotype value k. The 
females with k = 16 mate only with males of phenotype 
k = 16. If Y = 0.5, the female population is not selective 
in mating; that is, panmictic mating is the behaviour of 
the population. For Y = 1.0, the female population has 
a completely assortative mating behaviour. 



III. RESULTS 

We present now simulation results for which the val- 
ues of the parameters related to the first pairs of bitstring 
of the individuals' genomes were chosen to be: A = 10, 
6 = 5 and M = 1. The specialist populations have phe- 
notypes k G [0,ni = 13] and k G [n 2 = 19,32], Eq. 
(|2.ip and each female chooses among N m = 5 males. We 
start the simulations with all bitstrings randomly filled 



with zeroes and ones. The initial populations typically 
consist of 60000 individuals, half males and half females. 
The equilibrium population sizes depend on the carrying 
capacity, but are never smaller than 38000 individuals. 

A. Simulations for A. Labiatus species 

We took Z = in the resource distribution, F(j, t) of 
Eq. (|2.4[) . and chose for the mutation probability per lo- 
cus of the ecological and for the mating trait the values 
Alj = 1.0 and = 1.0, respectively. For all values of 
the sexual selection strength, Y, the phenotype frequency 
of the ecological trait is a stationary gaussian distribu- 
tion, Fig. [3(a). The phenotype frequency of the mating 
trait is a bimodal distribution for Y > 0.7 and a gaussian 
distribution, with a mean value k = 16, for the other Y 
values, Fig. [3(b). 

(b) Ecological character at X = 0.8, Z = and M k = 1.0 




(a) Mating character at X = 0.8, Z = and M k = 1.0 




FIG. 2: Results for the model without disruptive selection and 
with a competition strength X — 0.8. The phenotypes frequen- 
cies, for different sexual selection strengths, were measured 
during the last 10 6 simulation steps, (a) distribution of the 
phenotype, j, related to the ecological trait and (b) distribution 
of the phenotype, k, related to the mating trait. 

All these results, as well as the next one shown, are 
valid for all values of the competition strength smaller 
than one, X < 1. If we change N m — 5 to smaller values, 
the phenotype frequency related to the mating trait is no 
longer a stationary distribution. If the mutation proba- 
bility per locus of the ecological trait, Mj, is smaller than 



4 



1.0, for example Mj = 0.1, the distributions in Fig. [2] 
(a) and (b) do not change. When the mutation proba- 
bility of the mating trait, Mfe, changes from 1.0 to 0.1, 
for example, the phenotype frequency of the ecological 
trait, Fig. [2] (a), does not change either. However, for 
Mk = 0.1 and Y > 0.5 the phenotype frequency of the 
mating trait, Fig. [2](b), changes to a unimodal distribu- 
tion peaked at fc = or fc = 32, with a 0.5 probability 
for each. For example, the phenotype frequency of the 
mating trait, for Y > 0.5, is a distribution with a peak at 
k = 32, while for other values of Y and still Mk = 0.1 the 
distribution is bimodal peaking at k = 16 and k = 32, 
Fig. [3] That means, if the mutation probability of the 
mating trait, Mk, has small values, the population suffers 
a genetic drift which becomes faster as the value of Y be- 
comes larger, meaning that the behaviour of the female 
population is predominantly one of assortative mating. 
In other words, for strong assortative mating the disrup- 
tion of the mating trait is not favoured for small values 
of the mutation probability of this trait when disruptive 
natural selection, caused by resources distributions, is 
not present in the population. 



Mating character at X = 0.8, Z = 0.0 and M k = 0.1 
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B. Simulations for A. Zaliosus species 

The characteristics of the A. Zaliosus species is to have 
a unimodal distribution for the mating trait and a bi- 
modal distribution for the ecological one. We have al- 
ready seen, in the previous section, that for Z = there 
is no splitting of the ecological trait, unless if Y = 1.0 
and X = 1.0. On the other side, for Y — 1.0 and X = 1.0 
the mating trait is also splitted, which is not the case of 
A. Zaliosus species. So in order to simulate the A. Za- 
liosus process of speciation, we will take Z > and will 
first study each trait separately. 

(b) Ecological character at X = 0.5, Y = 0.5 and M k = 1.0 




(a) Mating character at X = 0.5, Y = 0.5 and M k = 1.0 




FIG. 3: The phenotypes frequencies of the mating trait for dif- 
ferent sexual selection strength, Y , without disruptive natural 
selection and at a small mutation probability of the mating 
trait, M k . 



FIG. 4: The phenotypes frequencies for different resources dis- 
tributions, Z , for Y = 0.5 and X — 0.5, (a) of the ecological 
trait and (b) of the mating trait. 



In what follows, the character determining sexual se- 
lection in the population is the colour of the individuals 
and the character that determines natural selection is the 
individuals' jaws morphology, which is the case in the Mi- 
das cichlid species complex. In common with A. Labiatus 
species, FigfTJ it is possible to find polychromatism and 
monomorphism in the ecological character for all values 
of the asymmetrical competition strength between the 
intermediate and specialist phenotypes smaller than one, 
X < 1 in Eq. (|23|) and (f2~2"|) , provided the following con- 
ditions are met: (i) no disruptive natural selection Z = 0, 
(ii) a sexual selection strength Y > 0.7 in the population, 
and (iii) a large value for the mutation probability of the 
mating character Mk = 1 , Fig. [J] (b) . 



Fig. H] (a) shows the distribution of the ecological trait 
for Y = 0.5 and X — 0.5. From this figure we see that 
this distribution changes from a unimodal one to a bi- 
modal distribution depending on the value of the dis- 
ruptive natural selection strength, Z. For Z » 1, the 
existing intermediate phenotypes belong to individuals 
that die before reaching the minimum reproductive age. 
Fig. [4] (b) shows the distribution of the mating trait, 
also for Y = 0.5 and X = 0.5. It can be seen that the 
distribution is unimodal only for Z < 1; for Z « 1 it 
is almost bimodal but not completely, since the interme- 
diate population, fc = 16, is appreciable. The existence 
of the intermediate phenotypes is due to the panmictic 
behaviour of the population, Y = 0.5 
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The phenotype frequency of the ecological trait, Fig. [4] 
(a), does not change if we vary the mutation probabilities 
0.1 < Mj < 1.0 and 0.1 < M k < 1.0. The phenotype 
frequency of the mating trait also does not change for 
0.1 < Mj < 1.0. The same is not true when < 1.0, 
since then the distribution is trimodal, Fig. [5l with peaks 
at k = 0, k = 16 and k = 32 when M k = 0.1. 



Mating character at X = 0.5, Y = 0.5 and M k = 0.1 
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FIG. 6: The phenotypes frequencies of the ecological trait in a 
case without disruptive natural selection and for strong sexual 
selection, Y = 1.0. 



IV. DISCUSSION 



FIG. 5: The phenotype frequencies of the mating trait for 
different resources distributions, Z , and at a small mutation 
probability of the mating trait, . 

From these results we may conclude that when Z ^ 0, 
the splitting of the mating trait is favoured by small mu- 
tation probabilities, Mk, of this trait, which is not what 
we want for the A. Zaliosus species. From Figs. |4] (a) 
and (b) we can see that the proper regions of parame- 
ters to simulate the speciation process of this species are 
Z > 0.4, bimodal distribution for the ecological trait, 
and Z < 0.8, unimodal one for the mating trait. 

C. Simulations for A. Citrinellus species 

One way to obtain the polycromatism and polymor- 
phism characteristic of the A. Citrinellus species, bi- 
modal distributions of the both traits, is to consider uni- 
form distribution of resources without disruptive natural 
selection. We have already shown that simulations with 
Z = 0.0 and Y = 1.0, only assortative mating, give a bi- 
modal distribution for the mating trait, Fig. H](b), inde- 
pendently of the competition strength, X, provided the 
mutation probability is = 1.0. However, to obtain 
also bimodal distribution of the ecological trait, with- 
out disruptive natural selection, it is necessary to have 
X = 1.0, that is, a symmetric competition between spe- 
cialist and intermediate phenotypes, as shown in Fig. [6l 

For Z > 0, it is also possible to split both phenotype 
distributions, but then it is necessary to consider small 
mutation probabilities of the mating trait, Mk, a weak 
assortative mating, Y > 0.5, and proper values of the 
resource distribution, Z, depending on the competition 
strength, X. For instance, Z > 0.7 for X = 0.5, FigH] 
(a). 



Although competition for resources and natural dis- 
ruptive selection appear together in Eq. (|2.2p and (|2.3[) . 
they lead to rather different situations. While competi- 
tion depends on the population sizes and affects equally 
individuals of the same phenotypic group, disruptive nat- 
ural selection caused by resources distributions acts on 
each particular individual, according to its phenotype. 
As a result, disruptive selection becomes more effective 
than competition. For instance, even for a panmictic be- 
haviour, disruptive natural selection may lead to a split- 
ting of both traits, depending only on the values of X and 
Z, as already pointed in A small mutation probabil- 
ity of the mating trait also favours this double splitting, 
as obtained in ]Vi\ . 

On the other hand, in order to split both phenotypic 
distributions without disruptive natural selection, Z = 0, 
it is imperative to have sexual selection, and now the 
splitting process depends on the values of Y and X, and 
in this case the mutation of the mating trait must be 
large in order to prevent genetic drift effects. 

Anyway, the A. Citrinellus case is the evolution of the 
splitting in both traits distributions must be driven by 
disruptive natural selection or by sexual selection, and 
the only statistical difference we found between these 
two scenarios is that the correlation between the traits is 
smaller when the splitting is driven by sexual selection. 
It so happens because the large mutation rate of the mat- 
ing trait, mentioned above, introduces large fluctuations 
in this correlation. 
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